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In an earlier work [S. A. Cheong and C. L. Henley, cond-mat/0206196 (2002)], we derived an exact 
formula for the many-body density matrix ps of a block of B sites cut out from an infinite chain 
of noninteracting spinless fermions, and found that the many-particle eigenvalues and eigenstates 
of pb can all be constructed out of the one-particle eigenvalues and one-particle eigenstates respec- 
tively. In this paper we improved upon this understanding, and developed a statistical-mechanical 
analogy between the density matrix eigenstates and the many-body states of a system of nonin- 
teracting fermions. Each density matrix eigenstate corresponds to a particular set of occupation 
of single-particle pseudo-energy levels, and the density matrix eigenstate with the largest weight, 
having the structure of a Fermi sea ground state, unambiguously defines pseudo-Fermi level. Based 
on this analogy, we outlined the main ideas behind an operator-based truncation of the density 
matrix eigenstates, where single-particle pseudo-energy levels far away from the pseudo-Fermi level 
are removed as degrees of freedom. We report numerical evidence for scaling behaviours in the 
single-particle pseudo-energy spectrum for different block sizes B and different filling fractions n. 
With the aid of these scaling relations, which tells us that the block size B plays the role of an inverse 
temperature in the statistical-mechanical description of the density matrix eigenstates and eigen- 
values, we looked into the performance of our operator-based truncation scheme in minimizing the 
discarded density matrix weight and the error in calculating the dispersion relation for elementary 
excitations. This performance was compared against that of the traditional density matrix-based 
truncation scheme, as well as against a operator-based plane wave truncation scheme, and found to 
be very satisfactory. 



I. INTRODUCTION 

An unsolved problem in many-body physics is to aug- 
ment numerical exact diagonalization, which is only fea- 
sible on comparatively small systems, so as to give the 
maximum information about the thermodynamic limit. 
We know from the renormalization group that only a few 
degrees of freedom, viz. those with low energies and long 
wavelengths, really matter. Hence it should be possible 
in principle to discard the irrelevant parts of the Hilbert 
space, but no method has been developed in real space 
for higher-dimensional, interacting lattice models. The 
density matrix renormalization group (DMRG)^"'* is the 
only successful method to date, but it is limited to one- 
dimensional chains, or two-dimensional systems that can 
be forced into one dimension as strips. ^'^ 

To get more mileage out of density matrix-based renor- 
malization group methods, surely we must develop a deep 
understanding of the structure of density matrices of the 
simplest possible systems, for which analytic results are 
available to guide us. Therefore, it is appropriate to be- 
gin by considering the ground state of a one-dimensional 
chain of noninteracting spinless fermions described by a 
nearest-neighbour hopping Hamiltonian 



Fermi sea ground state for a variety of translationally- 
invariant Hamiltonians with hoppings to further neigh- 
bors, provided their dispersion relation is monotonic so 
that the Fermi surface occurs at the same wavevector 
for a given filling. A 'block' is then identified within this 
one-dimensional system by choosing B sites that need not 
be contiguous, following which we can define the many- 
body density matrix ps of the block starting from the 
zero-temperature many-body wavefunction and tracing 
out all sites outside the chosen block of B sites. The 
basis of this paper is the simple factorized nature of ps 
and its eigenfunctions, which we derived exactly for non- 
interacting (spinless or spinfull) fermions/ by extending 
a technique introduced by Chung and Peschel.^ 

In Ref.8, it was shown that many-particle density ma- 
trix eigenstates are built from a set of single-particle cre- 
ation operators and eigenvalues, quite analogous to the 
energy eigenstates of a noninteracting system of fermions. 
For a block of B sites cut out from a larger overall sys- 
tem, there are B such pseudo- creation operators f^, and 
the block's Hilbert space is therefore spanned by all 2^ 
products of the //'s. The density matrix of the block can 
be written as 



H = - 



1 + 



(1.1) 



Pb =^"^exp(-i?), 



(1.2) 



Here we observe that, in one dimension, one has the same 



which is the exponential of a pseudo- Hamiltonian given 
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by 

B 

H = Y.^ifUu (1-3) 
1=1 

with single-particle pseudo-energy spectrum ipi. 

This paper is devoted entirely to noninteracting 
ferrnions, because tlie analytic results of Ref.7 permit 
many calculations that could not be done by numerical 
brute force in an interacting case, or can be carried out 
only by quite different methods (such as Monte Carlo). 
Furthermore, our hypothetical renormalization or pro- 
jection algorithm using DM truncation would first be 
applied to a well-understood system in which the low- 
energy excitations behave as noninteracting quasiparti- 
cles (as in a Fermi liquid or a d-wave superconducting 
phase). Our imagined numerical method in such a sys- 
tem would extract renormalized pseudo-creation opera- 
tors which are related to the single-particle opera- 
tors {//} in the same way that quasiparticle creation 
operators are related to bare fermion creation opera- 
tors. The pseudo-dispersion relation for the renormalized 
{//} is thus expected to scale in the same fashion as the 
pseudo-dispersion (pi encoded in (1.3) of the noninteract- 
ing fermions considered in this paper. 

Through the calculations in this paper we aim to un- 
derstand the analytic structure oi ps, and to begin to un- 
derstand the quantitative errors due to truncation. This 
includes the question: what is the proper measure of er- 
ror? The most familiar measure, the retained fraction 
of total density matrix weight, does not seem to be the 
best measure of error, as evidenced by the small errors 
obtained in Section VI for the calculations of the dis- 
persion relation. Another question to be investigated is 
whether — when we are severely truncating the Hilbert 
space, and attempting only to obtain the low-energy exci- 
tations the density matrix eigenstates are the optimal 
basis. As shown by the comparison in Section VII, they 
are certainly superior to the other plausible basis (plane 
waves). 

In Section II A, we summarize; first our niuc;h improved 
understanding of the analytic structure of the density 
matrix for noninteracting fermions (and by implication 
for any Fermi liquid) following from (1.2) and (1-3), and 
elucidate the statistical-mechanical analogy between the 
density matrix eigenstates and the many-body states of 
a system of noninteracting spinless fermions. In Sec- 
tion II B, we summarize the important results which we 
obtained in Ref.7, giving exact relationships between the 
block density matrix ps, the Green function matrix G 
restricted to the block, as well as their eigenstates and 
eigenvalues. Then in Section II C, we develop the main 
ideas behind a operator-based density m,a,trix (DM) trun- 
cation scheme based on the statistical-mechanical anal- 
ogy described in Section II A. 

The relation between the single-particle density matrix 
eigenstates and single- particle energy eigenstates of a sys- 



tem of noninteracting spinless fermions also suggests how 
the distribution of single-particle pseudo-energies ipi are 
expected to scale with the block size B. Numerically, a 
scaling relationship between ^pi and B was found indeed, 
for the overall chain at various fillings n. Our analytical 
results from Ref.7 shed light on this eigenvalue scaling in 
two ways. Firstly, as in Ref.S, they allow numerical study 
of the density matrix for system sizes so large that they 
would be inaccessible to any other techniques. Secondly, 
the exact connection of the block density matrix ps to 
the block Green function matrix G gives hints about the 
eigenvalue distribution. All of these will be discussed 
further in Section III, and the implications of the scal- 
ing behaviour of the single-particle pseudo-energies are 
discussed in Section IV, where we derived the asymp- 
totic behaviour, in the limit of infinite block sizes, of the 
largest density matrix weight and the truncated weight 
VKf, which is the sum of weights of the density matrix 
eigenstates retained in the operator-based DM trunca- 
tion scheme. 

Compared to the traditional density matrix truncation 
scheme used in the DMRG, our operator-based density 
matrix truncation scheme gives for the same number of 
density matrix eigenstates retained a slightly larger dis- 
carded weight e = 1 — Wt (see FIG. 8 in Section IV D). 
This quantity gives a 0(e) estimate as an upper bound 
— a worst case scenario — for the error incurred when 
computing the expectation of a most general observable. 
As with examples in numerical integration, the perfor- 
mance of an algorithm in integrating some classes of func- 
tions may be much better than that expected from the 
straightforward error analysis. In any case, we are not 
really interested in arbitrary observables, but rather, in 
the dispersion relation of elementary excitations, which 
we calculate in Section VI. The results are highly encour- 
aging: the dispersion relation calculated in the operator- 
based DM truncation scheme differ from the true disper- 
sion relation by an amount much smaller than what is 
suggested by the discarded weight. 

Besides quasiparticle dispersion relations, real space 
correlation functions are also interesting quantities to 
calculate, and these invariably depend on the real space 
structure of the many-body ground state wavefunction. 
Since this wavefunction is to be written in terms of the 
one-particle density matrix eigenfunctions, it is impor- 
tant to understand the real space structure of these as 
well. The one-particle density matrix eigenfunctions kept 
in our operator-based DM truncation scheme have spa- 
tial structures that are very similar to each other. In 
Section V we look into a representative one-particle den- 
sity matrix eigenfunction, the pseudo-Fermi eigenfunc- 
tion, for each B, and found that they also obey a univer- 
sal scaling relation. Then in Section VIC, we check how 
well such a truncated basis of one-particle density matrix 
eigenfunctions can approximate the true single-particle 
wavefunction at the Fermi level, which is a plane wave 
with wavevector kp- We find the approximation to be 
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good even when less than one quarter of the onc-particle 
density matrix eigenfunctions are kept. This is impres- 
sive, considering the fact that in the operator-based DM 
truncation scheme, the number of many-particle density 
matrix eigenstates thus represented by the one-particle 
density matrix eigenstates retained constitutes a minis- 
cule fraction of the total number of density matrix eigen- 
states. 

Finally, for systems where we know that the true 
single-particle wavefunctions are plane waves, it seems a 
priori plausible that a plane wave-based operator-based 
truncation scheme might outperform the operator-based 
density matrix truncation scheme. We look into this pos- 
sibility in Section VII, and find that while there are a few 
aspects in calculating the dispersion relation where the 
operator-has ed plane wave (PW) truncation scheme out- 
performs the operator-based DM truncation scheme, the 
overall performance of the PW scheme is inferior to the 
DM scheme. We then conclude in Section VIII by sum- 
marizing the important results obtained in this paper, 
and discuss where all these might fit into the numerical 
analysis of an interacting system. 



II. OPERATOR-BASED DENSITY MATRIX 
TRUNCATION 

A. Structure of Density Matrix Eigenvalues and 
Eigenstates 

Because the Hamiltonian in (1.1) conserves particle 
number, the eigenstates of pB have definite particle num- 
ber, and may be grouped into various P-particle sectors, 
where P = 0,l,...,i?. A consequence of our fundamen- 
tal formulas (1.2) and (1.3) is that every eigenstate of ps 
has the form 



written as 



\xl) = n /i io) ■ 



(2.1) 



Each eigenstate is specified by a list of pseudo-occupation 
numbers {ni}, where ni^ = 1 for the factors contained 
in (2.1), and is zero otherwise. Furthermore, the corre- 
sponding eigenvalue, the density matrix weight, is simply 
given by 



wl = S ^ exp 



E 

1=1 



where the quantity 



E 



(2.2) 



(2.3) 



appearing in the exponent is the total pseudo-energy. In 

terms of the single-particle pseudo-energies {(pi}, the nor- 
malization constant of the density matrix in (1.2) can be 



^ = ^ exp 



E 



(2.4) 



where the summation is over all 2^ combinations of oc- 
cupations. 

It is immediately clear from (2.2) that the density ma- 
trix eigenstate of maximum weight corresponds to the 
minimum total pseudo-energy. This is obtained by set- 
ting ni = 1 for ipi < and ni = for ipi > 0. In complete 
analogy to the real energy of a noninteracting system 
of fermions, we simply fill up the single-particle pseudo- 
energy levels (PELs) from the lowest up to a pseudo- 
Fermi level. The maximiim density matrix weight al- 
ways turns out to occur with the block fractional filling 
that is closest to the bulk filling of the Fermi sea ground 
state. More generally, the maximum-weight state in the 
P-particle sector is obtained by filling the states with 
the P lowest single-particle pseudo-energies. Finally, it 
is clear that the next-highest weights, or equivalently the 
next-lowest total pseudo-energies in the P-particle sector, 
are obtained by making particle-hole excitations involv- 
ing the PELs near to the last one filled. 

The above analogy may be extended to note that (1.2) 
is exactly the density matrix that would be obtained (see 
for example, Ref.9) at temperature T = lii H were the 
Hamiltonian. The reciprocal of the normalization con- 
stant ^ of Pb in (2.4) just corresponds to the grand par- 
tition function for the block of B sites. Among other 
things, this implies that (n;), the average particle num- 
ber in a particular PEL, has the functional form of the 
Fermi-Dirac distribution. 

We will actually apply this idea in a slightly different 
way, so as to relate the single-particle pseudo- energies tpi 
for different block sizes. If we were dealing with an ac- 
tual Hamiltonian, the dispersion relation would imply a 
density of states which would be multiplied by the sys- 
tem size to obtain the actual distribution of states. Our 
numerical scaling results in Section HI A confirm that 
pseudo-energies behave similarly to real energies. 



Relation of Pseudo-Energies to Block Green 
Function Matrix 



In our earlier work^ we obtain an exact formula 

H = -Y. [loge G{t - G)-\. c\c., (2.5) 



which, with (1.2), relates pB to the block Green function 
matrix G, whose matrix elements are = {c\c^) with 
i and j restricted to sites within the block. Clearly (2.5) 
becomes (1.2) when the pseudo-Hamiltonian is diagonal- 
ized. Also, (2.5) tells us that the quadratic form of H 
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in (2.5) and G are simultaneously diagonalizable. If we 
denote by 

lx;) = //|0) = ExKi)s-|0)' (2-6) 

the single-particle eigenstate of H with eigenvalue y;, 

then {xi(j)} is the eigenvector of G with eigenvalue 
A; . The single-particle pseudo-energies are related to the 
eigenvalues of G by 



loge 



A; 



1-A;^ 



or equivalently, 



A/ = 



exp ipi + l' 



(2.7) 



(2. 



i.e. the eigenvalues of G are the average pseudo- 
occupation numbers (n;). Note that we sometimes write 
ifi If {I, B) to make explicit the dependence on block 
size B. We will assume that tpi are ordered from the most 
negative to the most positive values. 

Another notable result that was derived in Ref.7 is that 



=2-1 = det(l -G)= JJ(1 - A;). 



(2.9) 



Along with (2.5) and (2.7), (2.9) comprises the final in- 
gredients that allow numerical computation of the den- 
sity matrix even in very large blocks. Aside from the 
possibilities of truncation, (2.1) through (2.9) have com- 
pletely reduced a 2^ x 2^ diagonalization problem into 
a B X B problem, a computational shortcut which allows 
numerical studies of large blocks. 



C. Recipe for Operator-Based Truncation 

The analytical structure of (1.3) hints at the proper de- 
sign of a truncation scheme. The retained Hilbert space 
of a block should not be the span of those density matrix 
eigenstates whose weight exceeds a cutoff. Instead, we 
should implement a 'consistent' truncation, such that the 
truncated Hilbert space consists of exactly 2'"°»'' states, 
built from all combinations of Imax pseudo-creation op- 
erators fl, . . . , , acting on a block 'vacuum state' 
|0)^, and satisfying fermion anticommutation relations. 

In the traditional density matrix-based truncation 
scheme used in DMRG, the recipe for truncation is to sort 
all density matrix weights in descending order, and then 
retain only the eigenstates associated with the weights 
that exceed a certain cut off. Let us refer to this as 
the weight-ranked DM truncation scheme. In light of our 
understanding of the structure of the many-body den- 
sity matrix presented in this paper, we can see that the 



weight-ranked DM truncation scheme will certainly re- 
tain the eigenstate with maximum weight, the pseudo- 
Fermi sea described in Section II A, along with eigen- 
states that are 'particle excitations', 'hole excitations' 
and 'particle-hole excitations' from the pseudo-Fermi sea. 
If we arrange the entire collection of many-particle den- 
sity matrix eigenstates into a state graph, then the state 
graph looks like a _B-dimcnsional hypercubic lattice near 
the pseudo- Fermi sea. What the weight-ranked DM trun- 
cation does in this state graph picture is to remove nodes, 
and in effect cut bonds, out from this hypercubic lat- 
tice, producing a subgraph that is much less connected 
and containing tenuous links. Because of this, when 
the Hamiltonian is projected onto the weight-ranked DM 
truncated basis, spurious interactions arc introduced. 

We can apply the pseudo-energy analogy in choosing 
how to truncate, given the form of the density matrix. It 
is familiar, in the tnmcation used in Fourier-spacc-bascd 
quantum renormalization groups, to discard all single- 
particle degrees of freedom except for a shell around the 
Fermi surface. In the same way, let us discard all op- 
erators fi as degrees of freedom, except those for which 
\ipi\ is less than some threshold <p*. For all other we 
'freeze' n; by setting n; to its ground state value 



ni 



1, for ipi < -Lf^; 
0, for ipi > ip^. 



(2.10) 



This choice gives the maximum density matrix weight, 
among the eigenstates having any particular set of n; 
for the retained single-particle pseudo-energy levels. The 
spirit of this truncation scheme is similar to that used 
in quantum chemistry, ^^^^^ except that the notion of 
a Fermi surface is more fuzzy in atoms and molecules. 
The idea that truncation consists of decreasing the thick- 
ness of a shell of wavevectors around the Fermi sur- 
face appeared in the original renormalization group for 
a quantum-mechanical solid-state problem (Anderson's 
poor man's RG for the Kondo problem). This obvious 
notion — that the action is around the Fermi surface 
— necessarily appears in every effective form of trunca- 
tion intended for a metal (see for example, Ref.l4 and 
Ref.l5, among others). However, to our knowledge all 
such schemes used a basis of plane waves or of energy 
eigenstates. Our variation uses PELs in analogy to the 
use of energy eigenstates in these previous problems. De- 
riving from the density matrix, it makes sense only in 
procedures that involving cutting a real-space block out 
of a larger system. 

Within this operator-based density matrix (DM) trun- 
cation scheme, we can define an effective Hamiltonian for 
the truncated Hilbert space, just by taking the matrix el- 
ements of the true Hamiltonian between all the retained 
states. Using an operator-based truncation, this will have 
a particularly clean form: first replace each creation op- 
erator c| by the equivalent combination of all {//}; then 

replace ///j 1 for each single-particle pseudo-energy 
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Lpi < —99, (these are frozen to be always occupied), and 
otherwise remove all terms involving the operators that 
are truncated. Thus, if the original Hamiltonian has at 
most Zjiiax-fermion terms, the same will be true for the 
truncated Hamiltonian. This prescription shows that 
such a truncation is possible for general models, once one 
knows the appropriate density matrix, biit in this paper 
we have applied it only to noninteracting models. 



III. SCALING BEHAVIOUR OF EIGENVALUE 
DISTRIBUTION 

Since the many-particle density matrix eigenvalues are 

built, according to (2.2), from the single-particle pseudo- 
energies, the latter are the focus of our numerical inves- 
tigations. Now if our entire system consisted of the block 
in a pure state at T = 0, then every eigenvalue A; of 
the block Green function matrix G, being the average 
pseudo-occupation number {ni) of a PEL, would either 
be one or zero. At T > 0, A; follow a Fermi-Dirac distri- 
bution. We will see later that cutting out a finite block 
from a T = system, by tracing over the environment of 
the block, has a similar effect on the eigenvalues of G as 
would taking T > when the block is the whole system. 

In a translationally invariant system with filling n (at 
T = 0), a fraction n of the eigenvalues of G arc one, while 
the rest are zero. Cutting out a block of length B must 
smooth out this step (much as having T > makes it into 
a Fermi-Dirac distribution), and we expect the transition 
from one to zero to occur over a fraction This 
guess was inspired by the analogy of pseudo-energy ipi 
in (1.3) to the real energy, which near the Fermi level 
scales linearly with wavevector ~1/B. This scaling 
suggests the conjecture of a scaling form for the single- 
particle pseudo-energy like ipi « Bf{l/B), and indeed we 
find below just such a scaling form. 



A. Pseudo-Energies and Pseudo-Occupation 
Numbers 

In this subsection we calculate numerically the eigen- 
values A; of the block Green function matrix G, and use 
(2.7) to compute the single-particle pseudo-energies ipi. 
For a chain of free spinless fermions in its ground state, 
the matrix elements of the block Green function matrix 
G are 



Gij — 



sm.'Kn\i — j\ 
ttH - i\ 



(3.1) 



where n is the filling fraction. FIG. 1 shows how Aj, 
the eigenvalues of G, are distributed for different filling 
fractions n and different block sizes B. 



For n 



our numerical studies suggest a scaling 



relationship of the form 

^{l,B) 




FIG. 1: Distribution of A; for different filling fractions and 
block sizes. In order to compare A; for different block sizes, 
the interval I € [1,-B] is rescaled such that {I — ^)/B £ (0, 1). 
With this rescaling, A = i always occur at [{I — \)/ B] = n. 



where x = [{I - ^)/B] - i, as shown in FIG. 2. There 
arc two points on FIG. 2 wc would like to note. First of 
all, with our choice of the scaling variable x, the scaling 
function f{x) always passes through the origin, i.e. 



/(O) = 0. 



(3.3) 



Secondly, from FIG. 2, we see that f{x) is an odd func- 
tion, i.e. 



f{-x) = -fix), 



(3.4) 



which is what we would expect from particle-hole sym- 
metry when the overall system is at half-filling, and f{x) 
has a finite positive slope at a; = 0, i.e. 

no) > 0. (3.5) 

Similar scaling behaviours of the form 

ipil,B,n)^Bf{n,x) (3.6) 

are found for all n, with the generic scaling variable 

x={l- If)/B, (3.7) 

where 



If = nB 



(3.8) 



plays the role of a Fermi wavevector, and a filling 
fraction-dependent scaling function f{n,x), as shown in 
FIG. 3. The scaling functions continue to satisfy 



and 



Bfix), 



(3.2) 



/(n,0) = 0, 

/'(n,0) >0, 



(3.9) 



(3.10) 
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FIG. 2: Plot of -B-MogJAi/(l - A()] as a function of the 
scaling variable x = [{I — — 1/2 for various block sizes at 
half-filling, showing a scaling collapse onto the scaling func- 
tion f{x). For B > 23, the largest and smallest pseudo- 
energies arc severely affected by numerical truncation errors 
in the diagonalization routines, and thus not shown. 



but f{n, x) is no longer an odd function of x for n ^ 
i. Instead, the particle-hole symmetry inherent in our 
model is manifested as 



f{n,-x) = -/(I - n,x). 
Prom (2.8) and (3.6), we can write A; as 

1 



Xi = 



exp(B/) + 1 ' 



(3.11) 



(3.12) 



which tells us that f{n,x) plays the role of a disper- 
sion relation €{k), while B plays the role of the inverse 
temperature /?. This confirms our suspicion that the ef- 
fect of cutting a block out of an overall system in its 
ground state at T = is to ascribe to the block an effec- 
tive temperature. As expected, this effective temperature 
approaches zero as the block size is increased, since we 
are keeping more and more information about the overall 
system, which we know to be at T = 0. 
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FIG. 3: Scaling collapses for (from top to bottom) n = 0.1, 
0.3, 0.7 and 0.9, plotted against the scaling variable x = {I — 
^)/B — n, for three different block sizes: B = 15 (o), B = 23 
(□) and B = 67 (o). 



Having understood that A; is related to the scaling 
function f{n,x) as in (3.12), we are now ready to in- 
vestigate the scaling behaviour of the normalization con- 
stant which is related to A/ by (2.9). We do this 
first at half-filling. As can be seen from FIG. 1, at 
half-filling roughly half of the A; are approximately one, 
whereas the other half are approximately zero. The prod- 
uct 11/(1 ~ '^i) is therefore determined primarily by the 
~B/2 Ai's that are nearly one. For these eigenvalues, 
cxp[Bf{x)] <C 1 and thus 1 — A; » exp[i3/(.T)] (when it 
is clear what the filling fraction is, we will drop the ar- 
gument n in f{n,x) to keep the notations neat). With 
this, we find that 

B B/2 

= 11(1 - AO « n ^MBfix)] 
1=1 1=1 

«exp|i?y°^ /(x)da;| (3.13) 

= exp|-By /(a;)rfa;|, 

where we made used the observed odd symmetry of the 
scaling function in (3.4) when the overall system is at 
half-filling, so that the integral within the exponent is 
positive. From (3.13), we see that decreases expo- 
nentially with block size B. 

In general, for n not too close to zero, where the argu- 
ment that those A;'s that are near one makes the domi- 
nant contribution holds, we find that 



exp \b f{n, x) dx 



exp 



j-B^ f{l-n,x)dxY 



(3.14) 
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where we made use of (3.11). For n very close to zero, 
there are a handful of Aj's of 0(1), and the rest are all 
nearly zero, behaving like A; w exp[— B/(n, x)]. Ignoring 
these handful of 0(1) Aj's, we find that the contribution 
to from those A; <C 1 is proportional to the product 
n,(l - AO « 1 - A;, and so 



oc 



^-Bf{n,x) 



(3.15) 



The integral can be evaluated as a cumulant expansion, 
but we can already see that for large B, the integral will 
not bo important, and thus derives most of its value 
from the few 0(1) A;'s. In contrast, when fi is very close 
to 1, then most of the eigenvalues A; of G are close to 1, 
and these continue to dominate the product ^((l — A;), 
and the asymptotic formula derived in (3.13) continues 
to be valid. 



A. The Gapless Chain of Noninteracting Spinless 
Fermions 

Instead of diving in to look at Wt or e, let us con- 
sider first a related question: how large is the maximum 
weight for a block of B sites embedded within an over- 
all system of gapless noninteracting spinless fermions? 
For our discussions, we will consider the half-filled case; 
it will be straightforward to extend the arguments pre- 
sented below to n ^ i. For convenience, let us take B 
to be even.^^ Let us denote by where F = B/2, 
the many-particle density matrix cigcnstatc having the 
largest weight. This state is always kept in the operator- 
based truncation scheme. Recall from Section II that, 
as shown in FIG. 4(b), this is the analog among density 
matrix eigenstates of the Fermi sea ground state among 
energy eigenstates. In this B/2-particle state, the single- 
particle pseudo-energy is filled up to just before x = 0, 
which we shall call the pseudo-Fermi level. The many- 



IV. DENSITY MATRIX WEIGHTS: 
IMPLICATIONS OF EIGENVALUE SCALING 

With our understanding of the structure of the many- 
particle density matrix eigenvalues and eigenstates de- 
veloped in Section 11 A, and on the scaling behaviour of 
the single-particle found in Section III, we want to now 
address the question of how much of the Hilbert space 
we can truncate. Clearly, the answer to this question de- 
pends on what measure of error we intend to use as our 
criteria for judging how well the truncated Hilbert space 
describes the physics associated with the parent model. 
In this section we look at the most common measure 
of error, used in the DMRG^'^ and quantum chemistry 
calculations:'^'^^^^ for a properly normalized density ma- 
trix, the density matrix weights wl satisfy the sum rule 



(4.1) 



If the ordinal numbers L are chosen such that wl is 
ranked in decreasing order, and a total of Lmax den- 
sity matrix eigenstates are retained, then the truncated 
weight 



Wt 



(4.2) 



L<L„ 



and the discarded weight 



e=l-Wt. 



(4.3) 



are frequently used as figures of merit for the truncation 
scheme, since for a bounded operator A, the truncation 
error in (A) is 0{e). 




(a) 



FIG. 4: Schematic diagram showing the three many-particle 

density matrix eigenstates, (a) |F — 1), (b) |_F) and (c) 
|F -|- 1), with the largest weights, for a block of B {B even) 
sites within a overall system that is half-filled. 

particle density matrix eigenstates of the block, of which 
there arc two, having the next largest weights will be 
called |F - 1) (FIG. 4(a)) and \F + 1) (FIG. 4(c)), hav- 
ing respectively one less or one more particle. 

We can understand the weight ratio wp-i/wp or 
Wp+i/wp as follows. By (2.2), \F) is the state with 
ni = 1 for Z = 1, B/2 and n; = for i = 
B/2 + I, B. The state \F + 1) differs only in 

having 71^/2-1-1 = 1) while \F — 1) differs only in hav- 
ing nB/2 = 0. Near the pseudo- Fermi level, the scal- 
ing fimction has a slope of /'(O), while the spacing be- 
tween adjacent pseudo-energies on the rescaled l/B axis 
is l/B. Thus — ~ /'(O)- But when the 

actual filling is n = 5, we know by particle- hole sym- 
metry that <Pb/2+i = -'4>B/2, so « -/'(0)/2 and 
'Pb/2+1 ~ /'(0)/2. It follows from (2.2), that 



exp(-/'(0)). 



(4.4) 



For n 7^ i, (3.6) would tell us that these ratios are ap- 
proximately exp(— /'(n, 0)). 

This is quite different from what would happen when 

the 'block' contains half of the entire system, as con- 
sidered in the standard DMRG algorithm, or in Ref.8. 
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If the fraction (B/N) in the block approached one, the 
state \F) must become the Fermi sea ground state of the 
overall system, and consequently contains all the weight. 
If the block is merely a finite fraction of the system, we 
still expect a much larger ratio than (4.4). It would be 
interesting to check the behaviour of the ratios in (4.4) 
for the case B, N ^ oo with B/N = 1/2, but we have 
not investigated this. 

B. The Gapped Chain of Noninteracting Spinless 
Fermions 

For the purpose of understanding the pseudo-energy 
spectrum of non- interacting systems better, we also con- 
sidered the dimerized tight-binding Hamiltonian 

N 

^ = E [1 + (-1)''^] + ' (4-5) 

where the hopping integral t is modulated by the (— 1)-'(5 

term to produce an energy gap. Henceforth we choose 
the scale of energy to be such that t — 1. This system 
was solved analytically by Gebhard et al,^^ wherein the 
Hamiltonian can be written as 

H= J2 E{k) (al^ak,+ - al_ak,-) , (4.6) 



single-particle pseudo-energies was found for all 5, each 
governed by a different scaling function f(fi,d,x). The 
scaling collapse plot for ^ = 0.5 is shown in FIG. 6, com- 
pared to the scaling function /(n, S = 0,x) for the gapless 
case. 



1.0 
0.5 



-1.0 









/W O-O 8 = 0.0 




/y 0-0 8=0.1 




/ <^8 = 0.3 








<-< 5 = 0.7 


1 , 1 


V-^ 5 = 0.9 

1,1,1 



[a-l/2)/B]-l/2 



FIG. 5: Plot of -S-MogJAi/(l - A;)] as a function of the 
scaling variable x = (/ — i)/B— i, for a block of size B — 23, 
with different hoping modulations 5 = 0.1,0.3,0.5,0.7,0.9. 
Pseudo-energies for |.x| > 0.2 for 5 > are not shown because 
these are severely affected by numerical errors incurred in 
the numerical integration and diagonalization routines. The 
various sets of straight line segments are intended to guide 
the eye in visualizing the data. 



with E{k) = ^/e^{k) + A^(k), where g(fc) = -2cosfc and 
A{k) = 2S sink. In terms of A(fc) and e(fc), we can define 
an angle (l>i~ such that tan20fc = A{k)/€{k), and whose 
sine and cosine wc denote as ak = cos 4>k, Pk = sin (j)k ■ In 
terms of these, the operators ak,+ and ak,- for the upper 
and lower bands respectively are given by 

Ofc,_ =afeC^-hi/3fcCfe+^, ^^^^ 
afc,+ = -/3fcCj, +iafeC^_,_^. 

At half-filling, the lower band is completely filled while 
the upper band is completely empty, and the ground state 
can written simply as 

l*)= n <-|0). (4.8) 

lfel<f 

For this ground state, the two-point function is given by 

G,, = ^ r dke^^^^-^^'^^tZ^^i^, (4.9) 
i-f Vcos2 k + S"^ sin^ k 

using which we can construct the block Green function 

matrix G, and hence using (2.7) the pseudoenergies which 
correspond to the density matrix eigenvalues. For a fixed 
block size oiB = 23, the pseudo-energy spectra for differ- 
ent hoping modulation 5 is shown in FIG. 5, compared 
to that of the gapless case. Scaling behaviour of the 



2 
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[(/-l/2)/B]-l/2 

FIG. 6: Plot of -B-MogJAi/(l - A,)] as a function of the 
scaling variable a; = (/ — \)/B — 1, for hoping modulation 
5 = 0.5 and different block sizes B = 9,13,17,23. The 
pseudo-cnorgics for \x\ > 0.3 arc not shown as these are 
severely affected by numerical errors in the numerical integra- 
tion and diagonalization routines. Also shown, is an approx- 
imate dashed curve for the scaling function f(x) for n = i 
and 5 = 0, obtained from the data for B = 23. The straight 
line segments are intended to guide the eye in visualizing the 
data. 

Repeating our analysis for the three density ma- 
trix eigenstates with the largest weights, we find again 
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that the ratios of density matrix weights wp+i/wp and 
Wf-i/wf to be independent of block size B when the 
overall system is at half-filling. However, these ratios de- 
pend strongly on the hopping modulation 5. As we can 
see from FIG. 5, the slope of the scahng curve at x = 
is steeper for larger 5. This indicates that — everything 
else being equal for finite B — a smaller fraction of den- 
sity matrix states is needed to capture the same total 
weight if the system is gapped. 

We have not investigated the case B/N 1/2, as 
in the standard DMRG algorithm, but we naturally ex- 
pect the ratio to increase in a gapped system. Thus |F) 
would be a better approximation to the ground state 
in a gapped system than in a gapless system, which is 
known as an empirical fact in the DMRG context. Our 
approach, if extended to the case B/N > 0, would give 
an analytic justification for this common observation. 

C. Largest Density Matrix Weight 




B 



FIG. 7: Plot of the largest density matrix weight wf as a 

function of the block size B for B even. The solid curve shown 
is a fit of the form wf ~ wf.oo + Awf exp(— B/Bq). The best 
fit to this small data set is obtained by neglecting the data 
points for i? = 2 and B = 4, for which we get wf,oo = 0.41, 
Awf = 0.26 and So = 0.11. 



For even B blocks on a gapless chain of noninteract- 
ing spinless fermions described the Hamiltonian (1.1), 
the largest density matrix weight wp can be numerically 
computed reliably till S « 20, and its dependence on B 
is shown in FIG. 7. Also shown in FIG. 7 is a fit of the 
numerical data to 



wf{B) = wf,oo + Awf exp(-B/So), 



(4.10) 



where wp.oo, Aw^? and Bq arc curve- fitting parameters. 
Here the exponentially decaying term is merely chosen to 
produce a good curve fit — we believe the S-dependence 
may be more complex — but what is interesting is the 
fact that Wf tends to a constant, wf.oo, in the limit of 
B oo. We find that we can understand this in terms 
of the scaling formulas developed so far. 

From Section II we saw that the largest many-particle 
density matrix weight wf corresponds to the situation 
for which all PELs below the pseudo- Fermi level (pF = 
are occupied, and all those above are empty. This means 
that 



WF = n 



(4.11) 



Using the fact that J2 can be written explicitly as 

^ = [](l + e-'^') ' (4-12) 



we then find that 



-<f>i 



Wf 



11 i + e-v^i l-|-e-*'i -'--•■1 

1<If 1>If I 



(4.13) 

To evaluate wp, we evaluate first its logarithm, which 
-log^w;F = X]loge(l + e"l*"l)- (4.14) 



Here we make two approximations. Firstly, because of 
(3.6), we know that ipi oc i?, and so except for a hand- 
ful of single-particle pseudo-energies ^pi very near ipp, all 
the exponentials are very small numbers. Using the ap- 
proximation logg(l -|- a;) ~ a; for x <C 1, we write (4.14) 
as 



- logg wpK^e- 



(4.15) 



Secondly, we note that because of (3.6), single-particle 
pseudo-energies far away from i^p = Q will contribute 
negligibly to the above sum. For B sufficiently large, 
those single-particle pseudo-energies making significant 
contribution in (4.14) will lie within a small interval 
about Ip where a linear approximation of the form 

^i^Bf'{0)^-^=f'mi-lF) (4.16) 

adequately describes the pseudo-dispersion relation. 
Substituting (4.16) into (4.15), we find then that 



B 



- logg ~ ^ 



e 
1=1 

oo 
1>If 



f'mi-h 



(4.17) 



-f'mi-iF) 



This is a geometric series which we can readily sum to 
give 

2 



■ loge Wp 



1 



(4.18) 



exp(-/'(0))' 

i.e. the largest density matrix weight wp is found to ap- 
proach a constant value of 



Wp = exp 



1 - exp(-/'(0)) 



(4.19) 
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as B ^ oo. From FIG. 2 and FIG. 3, we see that 
/'(O) 5, and so this asymptotic value of wp is aproxi- 
mately 0.13. This is smaller than the wf,oo = 0.41 found 
numerically. 



D. Discarded Weight 

Now that we understand more about the scaling be- 
haviour of the largest density matrix weight wp, let us 
analyze the discarded weight incurred by the operator- 
based DM truncation scheme. We compute numeri- 
cally the discarded weight incurred by the operator-based 
DM truncation scheme and that incurred by the weight- 
ranked DM truncation scheme, and show them in FIG. 8 
as a function of the number of many-body states kept 
as a comparison. As we can see, the discarded weight 
incurred by the operator-based DM truncation scheme 
is larger compared to the weight-ranked DM truncation 
scheme, for the same number of many-particle eigenstates 
kept. This is expected, since the wciiglit-ranked DM trun- 
cation scheme is by definition the most efficient scheme 
in exhausting the sum rule given in (4.1). In spite of 
this seemingly poorer 'convergence' property, we believe 
that the operator-based truncation scheme has advan- 
tages that cannot be reproduced by the weight-ranked 
truncation scheme, to be argued in detail in Section VI. 



— B = 1, weight-ranked 

— S = 13, weight-ranked 

— S = 15, weight-ranked 
• B = 1, operator-based 
□ B = 13, operator-based 
^ B = 15, operator-based 




FIG. 8: Discarded weiglit e as a function of the number 
of states kept: weight-ranked (as done in the DMRG), and 
operator-based. 

Writing the total density matrix weight explicitly as 

w=^-^j{{i+ e-'^o n (1 + «"'^') n + ^"'^')' 

kept below above 

(4.20) 

where the subscript 'kept', 'below' and 'above' refer to 
PELs retained, approximated as always occupied and ap- 
proximated as always empty in the operator-based trun- 
cation scheme respectively. The truncated weight Wt cal- 



culated from the operator-based truncation scheme is 
Wt = a-^\{{l + e-'^^) J] e-"^' J] 1. (4.21) 

kept below above 

Since = 1, the ratio Wt/W = Wt can be written as 



Wt= Y\ — — TT — - — 

below above 



(4.22) 



below 



above 



This expression has a simple interpretation in terms 
of the pseudo-occupation numbers {A;}. Using (2.8), we 
find that we can write Wt as 



Wt 



n A. 11(1-^0' 



(4.23) 



below above 



i.e. the truncated weight Wt is given by the product of 

pseudo-occupation numbers A; of those PELs we insist 
are always occupied, together with the product of the 
single-hole pseudo-occupation numbers (1 — A;) of those 
PELs we insist are always empty. From FIG. 1 we see 
that A; changes fairly rapidly from A; < 1 to A; > 0, over 
a small range of PELs. Therefore, it appears that there is 
a fairly large range of /'s for which A; is very close to one 
or very close to zero. However, this does not mean that 
we should perform an operator-based truncation scheme 
keeping only the small number of PELs whose Aj's are 
significantly different from one or zero. This is because 
Wt is bounded from above by 

Wt< Yl A„,ax n (1 - ^-i") ^ (A*)''-^''', (4.24) 



below 



above 



where 7 is the fraction of PELs retained in the operator- 
based truncation scheme, and 



A* — max(Amax) 1 ~ Amln)- 



(4.25) 



Because the exponent is 0{B), this number can still be 
very small. 

This brings us to the question we posed in the begin- 
ning of this section: how much of the Hilbert space do we 
truncate? If Wt is the only criterion then we see that a 
compromise is necessary. For a small block, the number 
of PELs with A; significantly different from one or zero is 
a sizeable fraction of the total number of PELs, but this 
number is manageable. For a large block, the number of 
PELs with A; significantly different from one or zero is 
a tiny fraction of the total number of PELs, but we still 
need 7 to be reasonably large for Wt to be appreciable in 
magnitude. This of course means that an unmanageably 
large number of PELs has to be retained. 

To make the above discussions more water-tight, let 
us make use of the scaling relations obtained thus far to 
find a formula relating the truncated weight Wt to both 
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the block size B and the fraction 7 of PELs retained. 
Taking the logarithm of (4.22) we find, using the fact 
that exp(— <C 1 for Z far below Ip, and exp(— t^;) <C 1 
for I far above Ip, that 



■log,W^t= J2 loge(l + e'^0+ E loge(l + e" 

below 



above 



below 



E 

KIf 



above 



E 



l=lF--yB/2 



E^-^'- E 



-V! 



E 

l=lp-^B/2 



E 



(4.26) 



where $* is a constant. 

If B is large and 7 small, then the linear approximation 
(4.16) for ipi is valid, in which case the two sums in (4.26) 
are equal, and given by 



If+'iB/2 



If+iB/2 



E 



E 



-f'(0){l-lF) 



1 - exp(-7i3/'(0)/2) 



(4.27) 



1 - exp(-/'(0)) 
With this, we can write Wj as 



Wt « W* exp 



? (l 

1 _ e-/'(0) 



-7B/'(0)/2^ 



, (4.28) 



where W* = cxp(— >!)*). We can find W* by taking the 
limit 7 ^ 0, in which case we retain no degree of free- 
dom in the PELs. Within the operator-based truncation 
scheme, this means that we insist all PELs below ipp to 
be always occupied and all those above ipp to he always 
empty, i.e. only the density matrix eigenstate with the 
largest weight is retained, and we should have 



Wt = WF = W*, 



(4.29) 



and so 
Wt ~ wp exp 



( 1 _ p-lBf'{0)/2\ 



. (4.30) 



This can be simplified further, using (4.19) to get 



Wt « exp 



-7S/'(0)/2 



(4.31) 



In the limit of 7 —> 1, we see from the above expression 
that Wt does not tend to one, but we understand that 
this is because the linear approximation (4.16) is only 
valid for a small range of PELs about ipp, i.e. only for 
small 7. In this regime, we may further approximate Wt 
as 



Wt ~ exp 



wp exp 



J (, _ f{ohB \ 

e-/'(0) V "2 J 

.no) , 



1 



-/'(o) 



(4.32) 



where 



In 



jB 



(4.33) 



is the number of PELs retained. As we can see, for small 
7, the truncated weight Wt increases exponentially with 
^max- Also, whenever (4.32) is valid, we get approxi- 
mately the same truncated weight Wt whether we use 
B = 100 and 7 = 0.2 or B = 200 and 7 = 0.1. We will see 
in Section VI that whenever the retained jB PELs lies 
within the regime where the pseudo-dispersion relation 
is linear, the truncation errors are essentially determined 
by Zmax = 7-B- 



V. SINGLE-PARTICLE DENSITY MATRIX 
EIGENFUNCTIONS 



A. A Priori Expectations 

As noted already in Section II C, in the many-body 
eigenstates with largest weights, all the very negative 
PELs will be occupied and all the very positive PELs 
will be empty. The only PELs with significant varying 
occupancy are those near the pseudo-Fermi level. 

By construction, the many-body density matrix eigen- 
states with large weights constitute the likely configu- 
rations of the block. The difference between the large- 
weight eigenstates of the P-particle and (P -|- l)-particle 
sectors of the density matrix is in the application of a 
creation operator such that the pseudo-energy \(pi\ is 
close to (fip = 0. In real space, it is likeliest that we can 
add a particle near the ends of the S-site block, for one 
can imagine that, in the first configuration, this particle 
was just past the end in an adjacent block, and we merely 
hopped it a short distance across the boundary to create 
the state of (P + 1) particles on the block in question. It 
follows that the single-particle eigenfunctions with single- 
particle pseudo-energies near the pseiido-Fermi level have 
their greatest amplitude near the block's boundaries. In 
other words, it is the sites near the end that are most 
correlated with the environmental information that we 
discarded. 
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B. General Features 

As noted earlier, the eigenstates of pB are all built up 
from the eigenstates |x;) of G, which are simultaneously 
the one-partiele eigenstates oi ps- As such, the effects 
of basis truncation, particularly in obtaining a truncated 
expansion of the target state j^*), must be understood 
in terms of the features of these one-particle eigenstates. 
The real-space features of \xi) can most easily be un- 
derstood in terms of the corresponding eigenfunctions 
Xi{j), where j = 1, . . . ,B are sites on the block. At half- 
filling, the probability densities |xKi)P exhibit particle- 
hole symmetry, as is shown in FIG. 9 for the example 
of -B = 9. In general, by node counting, we see that the 
sequence of B single-particle eigenfunctions are in one-to- 
one correspondence with the sequence of B plane wave 
states on the block, where the ordinal number I of the 
single-particle eigenfunctions is related to the wavevec- 
tor k of the plane wave states on the block by 




y 



FIG. 9: Probability density of the normalized one-particle 
eigenfunctions xi U) (plotted against the scaling variable y = 
— i ) /B) of pb on a block of B = 9 sites at half-filling, show- 
ing the particle-hole symmetry of the overall system. The 
subplots are arranged in order of increasing pseudo-energy. 



C. Scaling Behaviour 

For odd B, the pseudo-energy (p{B +i)/2 sits at the 
pseudo-Fermi level, and we may call the correspond- 
ing eigenfunction XfU) = X(B+i)/2(i) the pseudo-Fermi 
eigenfunction. The probability density associated with 
XfU) has nodes at every even j, as shown in in FIG. 10 
for the case of B = 23. The most prominent feature of 
the pseudo-Fermi eigenfunction, i.e. the amplitude be- 
ing strongest near the boundaries of the block, was first 



0.30 




y 

FIG. 10: Probability density function IxfO', B)P of the 
pseudo-Fermi eigenfunction for B = 23, plotted against the 
scaled variable (j — \)/B. 

observed by White. ^ This appears to be a generic fea- 
ture that occurs in both integrable and nonintegrable 
1-dimensional systems. Using the example of a chain 
of coupled harmonic oscillators, Gaite explained this 
"concentration of resolution of quantum states near the 
boundaries" as a simple consequence of angular quanti- 
zation of the density matrix. 

To analyze _B)p (where we write the B depen- 

dence of XF{j) more explicitly) more carefully, we first 
rescale the eigenvectors obtained from Octave^^ such that 

\Xf{{B+1)/2,B)\' = 1 (5.2) 

for B = 4p -Fl, p = 1, 2, . . . . For B = 4p -h 3, |xf((B + 
l)/2, B)]"^ = and the rescaling cannot be carried out as 
unambiguously as for the B = Ap+l series. This rescaling 
is harmless, since eigenvectors are only defined up to an 
arbitrary normalization. After this trivial rescaling, we 
find that the pseudo- Fermi probability density can be put 
in a scaling form 

\XF{j,B)f - N{B)giy)'s^^—p^, (5.3) 

sm wy 

where y = {j — ^)/B and g{y) is the scaling function 
shown in FIG. 11. 

In (5.3), N{B) is a i?-dependent normalization fac- 
tor, chosen to ensure that the pseudo-Fermi wavefunction 
given is properly normalized in the limit of B — > oo, i.e. 

B 

\imJ2\XF{j,B)f = l. (5.4) 

Although we cannot compute A''(B) analytically, we ven- 
ture a guess to its i?-dependence by noting that the func- 
tions g(%)) and sin Try are not very different from the func- 
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FIG. 11: Plot of the rcscalcd envelope function g(ij) for var- 
ious block sizes B = 4p + 1, p = 1,2, ... , compared against 
sin Try and 4y{l — y), where y = {j — ^)/B is the rescaled 
coordinate on the block. 



tion 4y(l — y), and so we expect 



B 



sin^ Try 



1 ^ 



j odd 



- + 1 
y i-y. 



[4y(i - 
1 



(5.5) 



which we can easily work out to have the form 



E 5 (y ) J ^^'^ ~ B (log, B + C) , 



sin^ ny 



(5.6) 



where C is a constant. Numerically, the best fit for N{B) 
in the range of block sizes B = 33 to B = 125 is obtained 
with 

N-\B) = 0.249BloggS + 0.668B. (5.7) 

Because of the enhanced amplitude near the edge of 
the block exhibited in the real-space structure of den- 
sity matrix eigenfunctions with single-particle pseudo- 
energies close to the pseudo-Fermi level, and conversely, 
enhanced amplitude near the center of the block exhib- 
ited in the real-space structure of density matrix eigen- 
functions with single-particle pseudo-energies far away 
from the pseudo-Fermi level, we worry that these eigen- 
functions might not be a good basis to use for expanding 
spatially uniform plane waves, which are the true single- 
particle energy eigenstates in our model. We address this 
concern in Section VI C. 



VI. OPERATOR-BASED DENSITY MATRIX 
TRUNCATION APPLIED TO DISPERSION 
RELATION 

In a gapless system, we conjecture that low-lying ex- 
citations above the ground state are built from the same 
operators as the long-wavelength fluctuations within the 
ground state. This supposition is certainly validated if 
the system has a continuous symmetry and the long- 
wavelength modes are Goldstone modes. In general it 
is justified by the relationships between correlation func- 
tions (for the ground-state fluctuations) and response 
functions (for low-energy excitations). 

Despite its poor convergence properties as far as ex- 
hausting the sum rule (4.1) is concerned, the operator- 
based truncation would still get the salient features of the 
physics right. We check this by projecting the Hamilto- 
nian in (1.1) onto the truncated set of fermion operators 
/(, and calculate the dispersion relation therefrom. There 
are two physical quantities of interest here: (a) for odd 
number of sites B, the middle band crosses the Fermi 
level, and we can ask how the Fermi velocity, given by 
the slope of the dispersion relation at the Fermi level, 
scales with B and the fraction 7 of fermion operators 
kept; or (b) for even B, a band gap develops as a result 
of truncation at the Fermi level, and we can ask how the 
size of this band gap depends on B and 7. 



A. Energy Gap at Fermi Level 



e(*), no truncation 
£(/:), 2 PELs truncated 
e(k), 4 PELs truncated 
E(*), 6 PELs truncated 




FIG. 12: Dispersion relation e(fc) for a block of _B = 8 sites, 
where the effects of truncating 2, 4 and 6 PELs are shown. 
For this block size, truncating 2, 4 and 6 PELs corresponds to 
fractions 7 = 0.75, 0.50 and 0.25 of PELs retained. For 7 = 
0.75, the energy bands (dashed curves) just below and above 
the Fermi level = agree with the true dispersion relation 
(solid curve) so well that the difference is not discernible at 
the scales presented in the figure. 

In FIG. 12, we show the general features of the disper- 
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sion relation e(fc) calculated within the operator-based 
truncation scheme, using the example of a block of i? = 8 
sites. Apart from the energy gap AiJ that opens up at the 
Fermi level e^;- = 0, we see that there is a one-to-one cor- 
respondence between the PEL truncated and the energy 
band absent from the dispersion relation. More precisely, 
if we order the energy bands and the PELs from the low- 
est to the highest as {ei(fe), . . . , esik)} and {ipi, . . . , tps}, 
then if we truncate PEL ipi, the energy band e;(fc) will 
also be removed from the numerically calculated disper- 
sion relation. For fixed 7, the gap AE decays exponen- 
tially with block size B, as is shown in FIG. 13, i.e. we 
have 



AE = AEo exp(-«;(7)B), 



(6.1) 



where K(-y) is an attenuation coefficient whose 7- 
depcndence is shown in FIG. 14. Here we sec also 
that AE{B,"f) for different 7 appears to converge onto 
AEo = AE{B = 0). Of course, there is no physical 
sense in talking about a block of zero size, but it is nev- 
ertheless a useful number to keep in mind when studying 
the scaling behaviour of AE{B,j) as 7 varies. AEq is 
approximately 4, which is the bandwidth of the exact 
dispersion relation, in all cases. 
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FIG. 13: Plot of the gap AE as a function of block size B for 
various constant fractions 7 = 1/4, 1/3, 1/2, 2/3, 3/4, 4/5 
of PELs retained. Also shown are fits to the data points for 
various 7 to AE{'y,B) = ASq exp(— /t(7)B). 

In particular, in the limit of 7 1, where all PELs 

are retained, the gap is exactly zero for all nonzero block 
sizes B. In this limit, if we start out at a 'gap' of AEq at 
a 'block size' of B = 0, then to have AE = at B = 1, 
we need the attenuation coefficient n to be infinite, i.e. we 
expect the limiting behaviour lim-y^i ^(7) = 00. On the 
other hand, in the limit of 7 ^ 0, where we retain none 
of the PELs, it is again physically meaningless to talk 
of a dispersion relation. Nevertheless, if we pretend that 
we are able to calculate a 'dispersion relation' in this 
limit, then it is reasonable, following the trend observed 
in FIG. 13, that the gap never closes, i.e. we expect the 
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FIG. 14: Plot of the attenuation coefficient ^(7) with which 
tlic spurious gap AE decays, as a function of the fraction 
7 of PELs retained. A cubic spline curve (dashed curve) is 
superimposed on the data points to aid visualization. Also 
shown (dotted curve) is — logg(l — 7). 



limiting behaviour lim..y^o '^(7) = 0- These limiting be- 
haviours appear to be borne out in the trend observed in 
FIG. 14. 

Another notable feature in FIG. 14 is the fact that 
^(7) « 7 for 7 <C 1, which is the regime we are most in- 
terested to apply the operator-based truncation scheme 
in. To appreciate the relevance of this observation, let 
us first note from FIG. 12 the general feature that the 
smaller the gap AE, the better the truncated disper- 
sion relation matches the true dispersion relation about 
the Fermi level. From (4.32) we saw that the truncated 
weight Wt depends only on the combination Zmax = 7-B 
in this regime, and as far as Wt is concerned, there is no 
difference whether we choose to keep 10 out of 100 PELs 
(7 = 0.1) or 10 out of 200 PELs (7 = 0.05). Here we see 
a similar exponential dependence on Zmax for the spuri- 
ous gap AE that arises due to truncation: if we write 
K « 7 in this regime, then AE « AEq exp(— B7) = 
AEq exp(— Zmax)- Such an exponential behaviour implies 
that we have very good control over the numerical ac- 
curacy of the dispersion relations — in particular near 
the Fermi level — calculated in the operator-based DM 
truncation scheme. 



B. Fermi Velocity 

When the block size B is odd, the central energy band 
crosses the Fermi level, and the quantity of interest be- 
comes the Fermi velocity vp- This can be determined 
from the truncated dispersion relation by taking the nu- 
merical central derivative of the central energy band at 
k = ±nTT/B. At half-filling, e(fc = ±7r/2B) = ex- 
actly because of particle-hole symmetry. This feature 
of the dispersion relation was found to be preserved in 
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the numerical dispersion relations computed within the 
operator-based DM truncation scheme. On the global 
scale, we find numerically that the shifts in the central 
energy band at the Brillouin zone center and Brillouin 
zone edge are such that vp > always. However, when 
B is large, the numerical diagonalization routine intro- 
duces artefacts on the energy scale of 10^^'^, resulting in 
the locally evaluated vp coming out to be very slightly 
less than 2. As such, we analyze the behaviour of vp as 
a function of B and 7 only for B < 150, as shown in 
FIG. 15. 




FIG. 15: Plot of Fermi velocity deviation {vf — 2) 

calculated from the truncated dispersion relation, as 
a function of the block size B, for various fractions 
7 = 1/15,1/13,1/11,1/9,1/7,1/5,3/15,3/13,3/11,3/9,3/7 
of PELs retained. Fits to average exponential decays are also 
shown. 

As can be seen from FIG. 15, the difference {vp — 2) 
decays more or less exponentially with B for various 7, 
i.e. 



{vp - 2) ^ exp(-^(7)B), 



(6.2) 



where ^(7) is the 7-dependent attenuation coefBcient for 
the average exponential decay. The 7-dependence of ^(7) 
is shown in FIG. 16. 



Real-Space Structure of Eigenfunctions of 
Truncated Hamiltonian 



The eigenfunctions of the untruncated Hamiltonian 
(1.1) are spatially uniform plane waves, with amplitude 
exp(ifcj)/\/S on site j of the block of B sites. These 
can be expanded in terms of the density matrix eigen- 
functions xiiji^)- Naively, we expect that if we drop 
those x.l{j^ B) associated with pscudo-cncrgics B) far 
from the pseudo- Fermi level ipp, as we would in our 
operator-based truncation scheme which removes these 
single-particle pseudo-energy levels as degrees of freedom, 
the remaining terms, all having enhanced amplitudes at 
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FIG. 16: Plot of the attenuation coefficient ^(7) as a function 
of the fraction 7 of PELs retained. Also shown (dashed line) 
is the expected behaviour ^ (7) = 7. 



the edge of the block, would sum to a function with en- 
hanced amplitude at the edge of the block. It would 

therefore scc!m like we are attempting to approximate a 
spatially-uniform plane wave with a function with the 
wrong real-space structure. 

Howcivcr, the key insight we gain from our study of 
block density matrices is that while the system- wide den- 
sity matrix po commutes with iJ in (1.1), the block den- 
sity matrix ps obtained by tracing down po does not 
commute with H{k), for all k. Therefore, after operator- 
based truncation H{k) — > H{k), we would need to di- 
agonalizc H{k) to find the truncated dispersion rela- 
tion ei{k). Thus, the function that would approximate 
the plane wave is not the latter's truncated expansion 
in terms of the eigenfunction of the one-particle block 
density matrix, but rather, a particular eigenfunction 
of H{k), which is an appropriate linear combination of 
the XiUjB) retained in the operator-based truncation 
scheme. We show in FIG. 17 the spatial structure of 
such a function, for various numbers of density matrix 
eigenfunctions retained. As we can see, for a block of 
B = 2'i sites, keeping 7 density matrix eigenfunctions 
with pseudo-energies around (fp would produce a decent 
approximation to the plane wave with k = 'k/2B. 



D. Discussion 

We speculate that the fact that the operator-based 

density matrix truncation scheme succeeds so well sug- 
gests that appropriate linear combinations of the density 
matrix eigenfunctions can closely approximate a plane 
wave with wavevcctor ±kp. This is only possible by tak- 
ing the difference of two eigenfunctions so as to cancel 
the enhancements of the envelope function seen at the 
block ends (see Section V). Indeed, the fact that we get 
the correct slope i;^ of the dispersion near kp suggests 
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FIG. 17: Amplitudes of eigenfunctions of the truncated 
Hamiltonian Hikw) for different numbers of PELs retained, 
for a block of B = 23 sites at half-filling. 



that by taking different weights, a continuously varying 
effective wavevector can be approximated. 

The fact that the goodness of approximation depends 
only on the number of eigenfunctions kept, means that 
wc approximate the wavefunction about as well in two 
successive blocks of B sites, as wc do in one big block 
of 2B sites. One could speculate that there might ex- 
ist some sort of approximate composition formula, anal- 
ogous to Clcbsh-Gordan formulas for combining angu- 
lar momenta, that provides the 2i?-site eigenfunctions in 
terms of the direct product of the B-site eigenfunctions. 



VII. OPERATOR-BASED PLANE WAVE 
TRUNCATION SCHEME 

As we saw in Section V, eigenstates of the density ma- 
trix pb are approximately plane waves (with wavevector 
Q determined by the boundary conditions on the block of 
B sites) modified by some envelope function. Apart from 
the effects of the envelope functions, the operator-based 
truncation scheme described above is likened to truncat- 
ing wavevectors Q far away from the Fermi wavevector 
kp . It is therefore natural to investigate how a operator- 
based truncation scheme based on plane waves would fare 
against that based on the density matrix eigenstates. 



A. Exact Dispersion at Zone Center 

Compared to the operator-has ed DM truncation 
scheme developed above, the most striking feature of 
the operator-based plane wave (PW) truncation scheme 
is that it gets the dispersion exactly right at the zone 
center, as shown for the case of B = 8 in FIG. 18, and 
for the case of B = 10 in FIG. 19. We understand this 
as follows: 




e(^), no truncation 
z(k), 2 PWs truncated 
e(fc), 6 PWs truncated 



FIG. 18: Dispersion relation e(fc) for a block of B = 8 sites, 
where the effect of truncating 2 and 6 jjlanc waves (PWs) 
are shown. For B = An, the Fermi level is located at the 
zone center, and e{k) is always gapless here regardless of the 
number of plane waves truncated. 



To evaluate the dispersion relation in a blocked de- 
scription, we start by defining the direct Bloch basis 
states 



1 



^NjB 



J 



\j,J), j = l,...,B, (7.1) 



where \ j, J) = c -_^jg |0) is the single-particle occupation 
number basis state at site j -t- JB along the chain. In 
this basis, the Hamiltonian (1.1) for a chain of N non- 
interacting spinless fcrmions take on a block diagonal 
form. Diagonalizing the B x B diagonal block 



H{k) 







-1 ••• -e-'^^J'^' 



-1 -1 
-10 



JkJB 







(7.2) 



for —tt/B < k < n/B then gives the dispersion relation 
within the reduced zone scheme. 

For the operator-based PW truncation scheme, we 
need to work with the plane wave states \Q,J) on each 
block of B sites, where the wavevector Q is determined by 
periodic boundary condition, i.e. exp{iQB) = 1. These 
plane wave states are related to the single-particle occu- 
pation number basis states by 



(7.3) 



A Bloch basis state parallel to (7.1) can be defined as 



\Q,k) 



1 



^e"=-^^|Q,J), (7.4) 
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where QB/2t: = 0, . . . , B - 1. Prom (7.3) and (7.4), it is 

easy to see that 

At the zone center, fc = 0, and the B x B Hamiltonian 
matrix in the | j, k) basis that we need to diagonahze be- 
comes 



H{Q) = 



0-10- 
-10-1- 
0-10- 

-10 0- 



-1 








(7.6) 



It is trivial to check that the eigenstates of this Hamil- 
tonian matrix are precisely the plane waves |Q, 0) on the 
block. Therefore, in the \Q,k) basis, H{k) is diagonal 
at A: = 0, and so truncating some plane waves from the 
Hilbert space produces no effect on the dispersion here. 

To be more precise, in performing truncation, a linear 
subspace of the Hilbert space is chosen, and the Hamil- 
tonian projected onto this subspace. If IV') is an eigen- 
statc of the full Hamiltonian, and if 1-0) is retained in the 
truncated Hilbert space, then it will continue to be an 
eigenstate of the truncated Hamiltonian, with the same 
eigenvalue. 



B. Energy Gap at Zone Boundary 

For even block sizes with B = An, the Fermi level is 
located at the zone center in the reduced zone scheme, 
and so there is no energy gap to speak of. On the other 
hand, for even block sizes with B = 4n+2, the Fermi level 
is located at the zone boundary. At the zone boundary, 
operator-based PW truncation introduces an energy gap 
AE at the Fermi level, as shown in FIG. 19 for B = 10. 

As in the case for the operator-based DM truncation 
scheme, we investigate the behaviour of the energy gap 
AE as a function of the block size B for a fixed frac- 
tion (1 — 7) of block states truncated. However, for the 
operator-based PW truncation scheme, the number of 
plane wave states that can be truncatcxi, if B is even, is 

4m -t-2,m = 0, 1,2, Thus the only realizable series 

of block sizes B on which we can perform fixed (1 — 7) 
truncation are of the form B = (j(4m + 2), g = 2, 3, 

The fraction 7 of block plane wave states retained 

is related to the series index q by 



7=1- 



(7.7) 



Half of these realizable series have block sizes that are 
multiples of 4, for which the Fermi level is at the Bril- 
louin zone center where the dispersion relation we have 




e(^), no truncation 
e(k), 2 PWs truncated 
elk), 6 PWs truncated 
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FIG. 19: Dispersion relation e(fe) for a block of B = 10 sites, 
where the effect of truncating 2 and 6 plane waves (PWs) 
are shown. For B = An + 2, the Fermi level is located at the 
zone boundary, and the operator-based plane wave truncation 
scheme introduces an energy gap AE here. 



shown in the previous subsection to be gapless. In this 
subsection we are interested in those block sizes for which 
q is an odd integer, since for these block sizes the Fermi 
level is at the Brillouin zone boundary, where a gap de- 
velops in the dispersion relation as a result of truncation. 
The behaviour of AE as a function of B for three series 
of q is shown in FIG. 20. 
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FIG. 20: Plot of AE as a function of block size B for g = 3, 
5 and 9, corresponding to the fractions 7 = 2/3,4/5,6/7 of 
block plane wave states retained. Also shown are the fits of 
the data points to the formula AE{B,^) = AEi{-y)/B. Prom 
the fits, wc have AEi = 2.58263 for 7 = 2/3, AEi = 1.57991 
for 7 = 4/5 and AEi = 1.13532 for 7 = 6/7. 



As can be seen from FIG. 20, the gap depends on block 
size as an inverse power law 



AE{B,j) 



B ' 



(7.E 
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where AEi (7) is a 7-dcpendcnt prefactor. This is in stark 
contrast to the exponential dependence (6.1) found for 
the case of the operator-based density matrix truncation 
scheme. 



C. Fermi Velocity 

For odd B, we again investigate the behaviour of vp 
as a function of B for the operator-based plane wave 
truncation scheme. The number of block plane waves 
that can be truncated is 4m + 3, m = 0, L 2, ... and the 
series of realizable block sizes are B = q{Arn+'i), g = 3, 5, 
.... Unlike in the operator-based DM truncation scheme, 
there appears to be two different systematic behaviours 
for 7), one for q ~ Ap—l and another for q = Ap+1 

{p = 1, 2, . . . ). We find that the Fermi velocity can be 
fitted very well to the formula 



Vf 



vf{i) + c+{i)/B, q = Ap-l; 
Vf{7) - C-{^)/B, q = Ap+l. 



(7.9) 
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FIG. 22: Plot of c± as a function of 7 in the operator-based 
PW truncation scheme for both the g = 4p — 1 series (c+(7)) 
and q = Ap + l (c-(7)) series. Both c+(7) and €-(7) appears 
to be converging towards 2. In fact, from the graph we find 
that |c±(7) - 2| fs 0.43(1 - 7). 



VIII. SUMMARY AND DISCUSSIONS 



The plots 01^^(7) and €±(7) as a function of 7 are shown 
in FIG. 21 and FIG. 22 respectively. 
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FIG. 21: Plot of as a function of 7 in the operator-based 
PW truncation scheme for both the q = 4p — l series and q = 
4p+l series. Also shown as the dashed curve is 2[1 — (1 — 7)^], 
which appears to fit the data points well near 7 = 1. 

As can be seen from FIG. 21, the exact value of the 
Fermi velocity is obtained only in the double limit of 
B — !■ 00 and 7^1. Compared to the operator-based 
DM truncation scheme, where we manage to achieve the 
exact Fermi velocity for any 7, this is clearly undesir- 
able. Furthermore, even very close to 7 = 1, the com- 
puted Fermi velocity approaches the limiting value vpi'y) 
as B~^. This is much slower than the exponential con- 
vergence of exp(— ^(7)5) found for the operator-based 
DM truncation scheme. 



A. Statistical-Mechanical Analogy and Scaling 

To summarize, we have in this paper developed an in- 
depth understanding of the structure of the eigenvalues 
and eigenstates of the density matrix pB of a block of B 
sites embedded within a infinite one- dimensional chain of 
noninteracting spinless fermions described by the set of 
fermion operators {c^,} with dispersion relation e(k). 

From Ref.7 we know that the block density matrix pB 
can be written in (1.2) as the exponential of a pseudo- 
Hamiltonian H given by (1.3) acting only within the 
block, which describes another system of noninteract- 
ing spinless fermions with fermion operators {fi}f^i and 
dispersion relation (pi. We use the prefix pseudo- when 
talking about operators {fi}f^i and energies fi to dis- 
tinguish them from the operators {c^,} and energies e{k) 
of the system we started out with. 

The many-particle eigenstates oi pB can then be in- 
terpreted simply as the many-body energy eigenstates 

(2.1) of the system of noninteracting spinless fermions de- 
scribed by H, and their associated density matrix weights 

(2.2) are then their statistical weights in the grand canon- 
ical ensemble, as though the block is at a finite temper- 
ature. Within this statistical-mechanical picture, we can 
apply intuitions learnt from the statistical mechanics of 
real fcrmionic systems, and talk about the filling of single- 
particle pseudo-energy levels (PELs) as dictated by the 
Fermi-Dirac distribution, the pseudo-Fermi sea with its 
pseudo- Fermi level tpp and particle- hole excitations of the 
pseudo-Fermi sea. 

Using results obtained in Ref.7, we make this 
statistical-mechanical picture more precise, by identi- 
fying the single-particle eigenstates \xi) in (2.6) of ps 
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as those built up from the eigenvectors {xi(j)}j j = 
1, . . . , B, of the Green function matrix G restricted to 
the block. The eigenvalue A; of G, related to the single- 
particle pseudo-energies ipi by (2.7) and (2.8), are the 
average pseudo-occupation numbers of the l^^ PEL. 

The statistical mechanics of real fermionic systems tells 
us that, at finite temperature, the physically important 
many-body states are those low-energy particle-hole exci- 
tations involving single-particle energy levels within ksT 
of the Fermi level. Single-particle energy levels far away 
from the Fermi level contribute negligibly to the thermo- 
dynamic properties of the fermionic system, and are pre- 
cisely the degrees of freedom to be truncated in a renor- 
malization group analysis. 

We capitalize on this insight, and described a recipe for 
operator-based truncation of the density matrix eigen- 
states, where out of the B pseudo-fermion operators 
{fi}E=i^ we truncate those /;'s associated with single- 
particle pseudo-energies (pi far away from ipp and retain 
^max of them with ipi « ipp- Within this operator-based 
truncation scheme, the effective pscudo-Hamiltonian act- 
ing on the truncated Hilbert space can be made to have 
the same form as the original pseudo-Hamiltonian, in the 
same spirit of renormalization group transformations in 
statistical mechanics or quantum field theory. 

Having laid out the basic principles behind our 
operator-based DM tnmcation scheme, we proceeded to 
look more closely into the distribution of single-particle 
pseudo-energies ipi, and how these scale with the block 
size B. There are two related questions that provide 
the motivation for doing this: (1) the statistical mechan- 
ics of real fermionic systems suggests that single-particle 
energy levels within /3~^ = ksT of the Fermi level are 
the physically relevant degrees of freedom — what then 
is the effective temperature (3^^ that we should use as 
the natural cutoff when performing operator-based trun- 
cation on (2) although we have associated the 
pseudo-dispersion relation ipi with the dispersion relation 
e(fc) of a real fermionic system of noninteracting spinless 
fermions, the wavevector k enumerating e(fc) is an inten- 
sive quantity whereas the ordinal number I enumerating 
ipi is an extensive quantity — what would the intensive 
analog of I that more closely parallels the wavevector k, 
and how would the pseudo-dispersion relation look like 
in terms of this intensive label? 

The natural answer to the second question would be to 
write the pseudo-dispersion relation ipi as a function of 
1/ B, totally analogous with how the wavevector k is enu- 
merated as 2-KmlN for a chain of N sites satisfying the 
Born-von Karman boundary condition. In fact, we find 
strong numerical evidence that suggests that the single- 
particle pseudo-energies, for various block sizes and fill- 
ings, satisfy a scaling relation of the form given in (3.6), 
where the scaling function f(fi,x) is the proper analog 
of the dispersion relation €{k), and the scaling variable x 
given in (3.7) is the proper analog of the wavevector k. 
Prom (3.12), we see that the block size B plays the role 



of inverse temperature. 

Our scaling results in Section 111 indicate that the den- 
sity matrix eigenstates and eigenvalues behave, as block 
size B is increased, very much as energy eigenstates and 
eigenvalues do when the system size is increased. In the 
latter case, we have a dispersion relation and are more or 
less sampling it at different wavcvcctors. It is not quite 
that simple in the density matrix case, since the scaling 
function (3.6) — our analog of the energy dispersion re- 
lation depends on the filling n. We only note that this 
analogy still lacks an analytical foundation. A more pen- 
etrating analysis is called for of the relation of G to ps 
(or equivalently, the effect on its eigenvalues of restricting 
G to sites on a block). 

Incidentally, we noted that our equation relating G to 
pB was valid at any temperature, but we assumed zero 
temperature throughout this paper. We expect nonzero 
temperature T would become a second scaling variable. 
Since T > has similar effects on the Green function 
G{r) as does the gap introduced in (4.5), we expect the 
scaling also behaves similarly and we did not investigate 
it. 



B. Scaling as a Guide to Truncation in Practice 

In a real application, it appears quite unlikely that 
Hamiltonians will be projected directly onto large blocks 
(meaning blocks of more than 16 sites). What then is 
the practical value of extracting scaling forms, if they 
are unambiguously seen only in blocks of 100 or more 
sites? One answer is that, even though it is an over- 
simplified cartoon for the non-asymptotic situations in 
which it usually gets applied, a scaling law is easier to 
grasp than brute numerical or graphical facts. 

The scaling relation (3.6) is a powerful tool that we 
can use to derive deeper understanding concerning the 
structure of the block density matrix, as well as various 
aspects of truncation. But in itself, the scaling relation 
provides only a partial answer to our first question, which 
is about how much of the Hilbert space of many-body 
states on the block of B sites we can truncate. To answer 
this question more completely, we looked at the three 
density matrix eigenstates \F — 1) and |F + 1) with 
the largest weights. Using the scaling relation (3.6), we 
find that the ratios wf±i/wf of their weights {wp being 
the largest density matrix weight) approach a constant, 
in (4.4) for a gapless system of noninteracting spinless 
fermions, as B — > oo. 

The same result was also found for a gapped system 
of noninteracting spinless fermions, for which we find 
scaling relations governed by gap-dependent scaling func- 
tions. Furthermore, the scaling relation (3.6) allowed us 
to conclude that as S ^ oo, the largest density matrix 
weight wp approaches a constant given in (4.19), and de- 
rive approximately the dependence on the number /max 
of PELs retained and the truncated weight Wt in the 
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operator-based DM truncation scheme. 



C. DMRG and Operator-Bsised Truncation 

It is difficult to coniparc our results with those ob- 
tained in the context of the DMRG, because that is an 
incremental method. Rather than obtain the density ma- 
trix of a large block all at once, each iteration of the 
DMRG takes an approximate density matrix for a block 
of B sites and produces an approximate density matrix 
for a block of {B + 1) sites. The fraction of weight kept, 
which is taken as the figure of merit, refers to the small 
truncation in each iteration. The cumulative DMRG 
truncation might be more appropriate to compare with 
our results for rather large blocks. 

Nevertheless, let us note that operator-based trunca- 
tion can be applied independent of whether we use an 
incremental or one-shot method. In particular, operator- 
based truncation could be used in a test run to ap- 
ply DMRG to a noninteracting Fermi chain. One is 
given a truncated list of t operators {/;} where / = 
{B -t+l)/2,...,{B + a-l)/2 for the original block, 
and a Hamiltonian projected onto it. One augments this 
list with the bare creation operators c'g^i on a new site 
that will be added, and defines the new Hamiltonian by 
adding the hopping to the new site. 

In light of the derivation of Ref.7, we anticipate that 
the density matrix of the augmented system's ground 
state must have the same quadratic form, with new oper- 
ators {//}, which could presumably be obtained merely 
by diagonalizing the single-particle sector. One simply 
deletes the least important member of this list to obtain 
a new truncated list, no longer than the original one. 

This difficulty notwithstanding, we still carried out a 
naive comparison of the performance of the operator- 
based DM truncation scheme against the traditional 
weight-ranked DM truncation scheme used in the 
DMRG, using the ability to exhaust the sum rule (4.1) 
for a given total numer Lmax of density matrix eigen- 
states retained as a criterion. The conclusion: while the 
operator-based DM truncation scheme do not exhaust 
the sum rule (4.1) as rapidly as the weight-ranked DM 
truncation scheme, Wt is still of 0(1), i.e. the significant 
parts of the total density matrix weight are 'captured' by 
the operator-base DM truncation scheme. 



1. Truncation and Dispersion Relations 

However, we believe it is more important to check how 

well a truncation scheme do by calculating physical quan- 
tities, rather than rely solely on the truncated weight Wt 
as a performance indicator. To this end, we calculated 
the dispersion relation of elementary excitations within 
the operator-based DM truncation scheme (an easy thing 



to do), and found that the error incurred decays expo- 
nentially as ^max, the number of PELs retained when 
^max B. This error is much smaller than 0(e), which 
is expected from a naive analysis based on the discarded 
weight e = 1 — Wt- Here, there is subtle worry that it 
may be that truncation works especially well for our cho- 
sen hopping Hamiltonian (1.1) is so local. A Hamiltonian 
with longer range hopping would have the same Fermi sea 
and hence the same density matrix, but the truncation 
errors might be worse. 

It would be desirable to also calculate the dispersion re- 
lation within the weight-ranked DM truncation scheme, 
and compare the results to those obtained within the 
operator-based DM truncation scheme. However, in the 
latter case it is problematic even to define the question, 
since each retained density matrix eigenstate would be 
a many-particle state. The new truncated Hamiltonian 
might be conveniently expressed in terms of the pseudo- 
creation operators {/;}, but many combinations of occu- 
pations would not exist. 

The situation would be somewhat analogous to taking 
a simple, noninteracting hopping Hamiltonian for spin- 
full fermions, and imposing a Gutzwiller projection (no 
doubly occupied sites). In effect, the projection made a 
noninteracting model into an interacting one. Similarly a 
weight-ranked truncation must introduce spurious inter- 
actions. Hence, even if a system containing several blocks 
were to be exactly diagonalized (using the truncated ba- 
sis) we could not immediately identify the elementary ex- 
citations. One would need, for example, to numerically 
compute a spectral function S{q^uj), where {hq,hw) are 
momentum and energy, and then locate peaks as a func- 
tion of huj. On the other hand, a system which is trun- 
cated using the operator-based truncation scheme can 
still be represented by a set of creation and annihilation 
operators. 

2. Inadequacy of Plane- Wave Truncation Scheme 

Following this, we argued, based on the real-space 
structure of the one-particle density matrix eigenfunc- 
tions shown in Section V, that an operator-based trunca- 
tion scheme can also be defined naturally using the basis 
formed by single-particle plane wave (PW) states on the 
block of B sites. The dispersion relation was calculated 
within this operator-based PW truncation scheme, and 
compared to that calculated within the operator-based 
DM truncation scheme. We find that, other than getting 
the dispersion exactly right at the zone center of the re- 
duced Brillouin zone, the operator-based PW truncation 
scheme is generally inferior to the operator-based DM 
tnmcation. Instead of decaying exponentially as Zmax, 
the error in the dispersion relation calculated within the 
operator-based PW truncation scheme decays as a power 
law which means that more single-particle basis 

states need to be retained in the operator-based PW 
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truncation scheme as compared to the operator-based 
DM truncation scheme. 



D. Towards Interacting Systems in Dimensions 
d>2 

It is not much harder in principle to analyze the pro- 
posed operator-based density matrix truncation scheme 
for noninteracting fermions in two dimensions. However, 
it will be harder to understand the scaling since we can- 
not simply rank the eigenvalues. One dimension was spe- 
cial because we know that each successive state has one 
more node than the previous one. A further very impor- 
tant difference is that in = 1 there arc just two Fermi 
points, whereas in d > 2 there is a Fermi surface. Thus, 
whereas in d = 1 the density matrix eigenstates near the 
pseudo-Fermi level arc related to the energy eigenstates 
at ikp, in d > 2 these eigenstates will correspond to 
mixtures of wavevectors from every piece of the Fermi 
surface. 

Obviously, noninteracting systems do not require nu- 
merical studes, so we must clarify how our results are 
relevant to the problem of interacting systems. Firstly, 
many (gapless) systems of interest are in a phase — Fermi 
liquid, d-wave superconductor — which are noninteract- 
ing in the low-energy, long- wavelength limit. When ap- 
plied to a Fermi liquid system, we expect (to the ex- 
tent that the truncation has separated out the low en- 
ergy modes) that any iterative renormalization scheme 
will converge on the noninteracting limit, and it must 
behave properly in that limit to have even the hope of 
success. Hence, for a density matrix-based scheme, the 
first order of business is to study the density matrix for a 
free fermion ground state (as in this paper) or for a BCS 
state, ^ i.e. that the density matrix will actually have the 
simple operator-based form, and hence its many-particle 
eigenstates of the density matrix really are built from the 
single-particle eigenstates. 

For an interacting system, the truncated basis should 
of course be constructed using the many-body density 
matrix for that system (not the noninteracting system). 
We expect that the operators generating this many- 
particle truncated basis will not just be those that cre- 



ate the 1-particle density matrix eigenstates (they were 
in the noninteracting case studied in this paper). More 
thought will be needed as to discover the best recipe to 
optimize the truncation rule so as to balance the needs 
of sectors with different particle numbers in a strongly 
interacting system. The simple algebraic structure of the 
noninteracting density matrix eigenstates, considered in 
this paper, has motivated investigations (in progress) of 
the relationships among the many-particle density matrix 
eigenstates for an interacting system. 

A separate reason why our results for the noninter- 
acting fermions are relevant to the study of interacting 
systems is that the scaling behaviour of tlw^ noninter- 
acting density matrix should be a good guide to that of 
interacting systems, although details may differ. This is 
in the same sense that mean-field theory is a good guide 
to the overall pattern of critical phenomena. 

However, other interacting models of interest sit at 
quantum critical points that are not described by quasi- 
particle interactions, or possess fractionalized excita- 
tions. Since we do not presently understand the proper 
way to write their wavefunctions in terms of a spatially 
blocked basis, nor the proper renormalization step to 
capture the interblock correlations in the fractionalized 
systems, we do not know if a plain block density ma- 
trix gives the proper basis for truncation of the states. 
Furthermore, in the absence of an analytic construction 
of the block density matrix, for example, for Laughlin's 
quantum Hall wavefunction or the one-dimensional Su- 
Schrieffer state, wc cannot proceed to scaling studies of 
large blocks like those found in the present paper, but 
numerical studies of such density matrices might be an 
illuminating subject for future research. 
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